################################################################################
# Joshua C. Fjelstul, Ph.D.
# "How the Chamber System at the CJEU Undermines the Consistency of the 
# Court's Application of EU Law"
# Journal of Law and Courts
# replication code for the computational formal model
################################################################################

# install ggminimal package (a ggplot theme)
# devtools::install_github("jfjelstul/ggminimal")

# load packages
library(tidyverse)
library(ggplot2)
library(ggminimal)
library(patchwork)

##################################################
# function to simulate a case
##################################################

simulate_case <- function(case_id, positions, facts, count_judgments, chamber_size) {
  
  # number of judges
  judges <- length(positions)
  
  # data (chamber assignment)
  panels <- matrix(data = 0, nrow = count_judgments, ncol = judges)
  for(judgment in 1:count_judgments) {
    panel <- sample(1:judges, size = chamber_size, replace = FALSE)
    panels[judgment, panel] <- 1
  }
  
  # simulate case facts
  mean_position <- NULL
  for(judgment in 1:count_judgments) {
    judge_indexes <- which(panels[judgment,] == 1)
    mean_position[judgment] <- mean(positions[judge_indexes])
  }
  
  # probability matrix
  probability_matrix <- matrix(data = NA, nrow = count_judgments, ncol = judges)
  for(judgment in 1:count_judgments) {
    for(judge in 1:judges) {
      
      # if the judge is on the panel...
      if(panels[judgment, judge] == 1) {
        
        # calculate the probability the judge votes yes
        probability_matrix[judgment, judge] <- pnorm(facts - positions[judge])
      }
    }
  }
  
  # outcome probabilities
  probability <- NULL
  for(judgment in 1:count_judgments) {
    chamber_probabilities <- probability_matrix[judgment,]
    chamber_probabilities <- chamber_probabilities[!is.na(chamber_probabilities)]
    probability[judgment] <- 1 - poibin::ppoibin(((chamber_size + 1) / 2) - 1, chamber_probabilities)
  }
  
  # create output
  output <- dplyr::tibble(
    case_id = case_id,
    judgment_id = 1:count_judgments,
    chamber_size = chamber_size,
    probability = probability,
    facts = facts,
    mean_position = mean_position
  )
  
  # return
  return(output)
}

##################################################
# function to run the simulation
##################################################

run_simulation <- function(count_judges, count_cases, count_judgments, mu_positions, sigma_positions, mu_facts, sigma_facts, chamber_sizes) {
  
  # create an empty list to hold output
  output <- list()
  
  # loop through cases
  for (i in 1:count_cases) {
    
    # message
    message <- stringr::str_c("\rdownloading ", i, " of ", count_cases)
    message <- stringr::str_pad(message, side = "right", pad = " ", width = 100)
    cat(message)
    
    # draw judge positions
    positions <- rnorm(count_judges, mu_positions, sigma_positions)
    
    # draw case facts
    facts <- rnorm(1, mu_facts, sigma_facts)
    
    # create an empty list
    output_case <- list()
    
    # loop through panel sizes
    for (j in 1:length(chamber_sizes)) {
     
      # simulate outcomes by panel size
      output_case[[j]] <- simulate_case(
        case_id = i,
        count_judgments = count_judgments,
        positions = positions,
        facts = facts,
        chamber_size = chamber_sizes[j]
      ) 
    }
    
    # stack tibbles for each chamber size
    output[[i]] <- output_case
  }
  
  # stack tibbles for each case
  output <- dplyr::bind_rows(output)
  
  # return Monte Carlo output
  return(output)
}

##################################################
# function to summarize the results
##################################################

summarize_results <- function(out_0) {
  
  # define pipe locally
  `%>%` <- dplyr::`%>%`
  
  # aggregate results by panel size and iteration
  out_0 <- out_0 %>%
    dplyr::group_by(chamber_size, case_id) %>%
    dplyr::summarize(variance = sd(probability) ^ 2) %>%
    dplyr::ungroup()
  
  # aggregate results by panel size
  out_0 <- out_0 %>%
    dplyr::group_by(chamber_size) %>%
    dplyr::summarize(mean_variance = round(mean(variance), 5)) %>%
    dplyr::ungroup()
  
  # return the results
  return(out_0)
}

##################################################
# function to calculate comparative statics
##################################################

calculate_cs <- function(sim_0, sim_1 = NULL, size_0, size_1, magnitude = FALSE) {
  
  # if one simulation is provided
  if(is.null(sim_1)) {
    
    # summarize the results
    res_0 <- summarize_results(sim_0)
    
    # get the variance for each condition
    variance_0 <- res_0$mean_variance[res_0$chamber_size == size_0]
    variance_1 <- res_0$mean_variance[res_0$chamber_size == size_1]
    
    # calculate the first derivative
    b <- variance_1 - variance_0
  }
  
  # if two simulations are provided
  else {
    
    # summarize the results
    res_0 <- summarize_results(sim_0)
    res_1 <- summarize_results(sim_1)
    
    # get the variance for each condition
    variance_0_0 <- res_0$mean_variance[res_0$chamber_size == size_0]
    variance_0_1 <- res_0$mean_variance[res_0$chamber_size == size_1]
    variance_1_0 <- res_1$mean_variance[res_1$chamber_size == size_0]
    variance_1_1 <- res_1$mean_variance[res_1$chamber_size == size_1]
    
    # calculate the first derivatives
    delta_variance_0 <- variance_0_1 - variance_0_0
    delta_variance_1 <- variance_1_1 - variance_1_0
    
    # calculate the second derivatives
    if (magnitude == FALSE) {
      b <- delta_variance_1 - delta_variance_0
    } else if(magnitude == TRUE) {
      b <- abs(delta_variance_1) - abs(delta_variance_0)
    }
  }
  
  # return the effect
  return(b)
}

##################################################
# run simulations
##################################################

# set seed
set.seed(12345)

# simulation 0: baseline
sim_baseline <- run_simulation(
  count_judges = 27, chamber_sizes = c(1, 3, 5, 7, 15, 27),
  count_cases = 1000, count_judgments = 1000,
  mu_positions = 0, sigma_positions = 1,
  mu_facts = 0, sigma_facts = 1
)

# simulation 1: judge heterogeneity
sim_judges <- run_simulation(
  count_judges = 27, chamber_sizes = c(1, 3, 5, 7, 15, 27),
  count_cases = 1000, count_judgments = 1000,
  mu_positions = 0, sigma_positions = 2,
  mu_facts = 0, sigma_facts = 1
)

# simulation 3: case facts
sim_case_facts <- run_simulation(
  count_judges = 27, chamber_sizes = c(1, 3, 5, 7, 15, 27),
  count_cases = 1000, count_judgments = 1000,
  mu_positions = 0, sigma_positions = 1,
  mu_facts = 1, sigma_facts = 1
)

##################################################
# save results
##################################################

# save
save(sim_baseline, sim_judges, sim_case_facts, file = "replication-materials/comparative_statics.RData")

# load
load("replication-materials/comparative_statics.RData")

##################################################
# summarize results
##################################################

res_baseline <- summarize_results(sim_baseline)
res_judges <- summarize_results(sim_judges)
res_case_facts <- summarize_results(sim_case_facts)

##################################################
# comparative statics
##################################################

calculate_cs(sim_0 = sim_baseline, size_0 = 27, size_1 = 15)
calculate_cs(sim_0 = sim_baseline, size_0 = 15, size_1 = 7)
calculate_cs(sim_0 = sim_baseline, size_0 = 7, size_1 = 5)
calculate_cs(sim_0 = sim_baseline, size_0 = 5, size_1 = 3)
calculate_cs(sim_0 = sim_baseline, size_0 = 3, size_1 = 1)

# positive effect

##################################################
# Figure 2, Panel A
##################################################

# plot data
plot_data <- tibble(
  x = res_baseline$chamber_size,
  y = res_baseline$mean_variance
)

# make plot
panel_a <- ggplot(plot_data) + 
  geom_segment(aes(x = res_baseline$chamber_size[1], y = -Inf, xend = res_baseline$chamber_size[1], yend = res_baseline$mean_variance[1]), linetype = "dashed", size = 0.25, color = "gray80") + 
  geom_segment(aes(x = res_baseline$chamber_size[2], y = -Inf, xend = res_baseline$chamber_size[2], yend = res_baseline$mean_variance[2]), linetype = "dashed", size = 0.25, color = "gray80") + 
  geom_segment(aes(x = res_baseline$chamber_size[3], y = -Inf, xend = res_baseline$chamber_size[3], yend = res_baseline$mean_variance[3]), linetype = "dashed", size = 0.25, color = "gray80") + 
  geom_segment(aes(x = res_baseline$chamber_size[4], y = -Inf, xend = res_baseline$chamber_size[4], yend = res_baseline$mean_variance[4]), linetype = "dashed", size = 0.25, color = "gray80") + 
  geom_segment(aes(x = res_baseline$chamber_size[5], y = -Inf, xend = res_baseline$chamber_size[5], yend = res_baseline$mean_variance[5]), linetype = "dashed", size = 0.25, color = "gray80") + 
  geom_segment(aes(x = res_baseline$chamber_size[6], y = -Inf, xend = res_baseline$chamber_size[6], yend = res_baseline$mean_variance[6]), linetype = "dashed", size = 0.25, color = "gray80") + 
  geom_point(aes(x = x, y = y), color = "black", size = 3) + 
  scale_x_continuous(breaks = c(1, 3, 5, 7, 15, 27)) + 
  theme_minimal() + 
  titles_minimal(title = "Panel A", x = "Chamber size", y = "Average variance")

##################################################
# Figure 2, Panel B
##################################################

# plot data
plot_data <- tibble(
  x = factor(
    x = c(1, 2, 3, 4, 5), 
    levels = c(1, 2, 3, 4, 5), 
    labels = c("27 to 15", "15 to 7", "7 to 5", "5 to 3", "3 to 1")
  ),
  y = c(
    calculate_cs(sim_0 = sim_baseline, size_0 = 27, size_1 = 15),
    calculate_cs(sim_0 = sim_baseline, size_0 = 15, size_1 = 7),
    calculate_cs(sim_0 = sim_baseline, size_0 = 7, size_1 = 5),
    calculate_cs(sim_0 = sim_baseline, size_0 = 5, size_1 = 3),
    calculate_cs(sim_0 = sim_baseline, size_0 = 3, size_1 = 1)
  )
)

# make plot
panel_b <- ggplot(plot_data) + 
  geom_point(aes(x = x, y = y), color = "black", size = 3) + 
  geom_segment(aes(x = x[1], y = -Inf, xend = x[1], yend = y[1]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[2], y = -Inf, xend = x[2], yend = y[2]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[3], y = -Inf, xend = x[3], yend = y[3]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[4], y = -Inf, xend = x[4], yend = y[4]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[5], y = -Inf, xend = x[5], yend = y[5]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_point(aes(x = x, y = y), color = "black", size = 3) +
  theme_minimal() + 
  titles_minimal(title = "Panel B", x = "Change in panel size", y = "Effect of change\non the average variance")

##################################################
# Figure 2
##################################################

# combine plots
figure_2 <- panel_a + panel_b + plot_layout(ncol = 2)

# save plot
ggsave(plot = figure_2, filename = "replication-materials/figure-2.pdf", device = "pdf", width = 10, height = 5, scale = 1.25)

##################################################
# heterogeneity
##################################################

# going from less heterogeneity to more heterogeneity
calculate_cs(sim_0 = sim_baseline, size_0 = 5, size_1 = 3) # cut-point standard deviation at 1
calculate_cs(sim_0 = sim_judges, size_0 = 5, size_1 = 3) # cut-point standard deviation at 2
calculate_cs(sim_0 = sim_baseline, sim_1 = sim_judges, size_0 = 5, size_1 = 3, magnitude = FALSE)
calculate_cs(sim_0 = sim_baseline, sim_1 = sim_judges, size_0 = 5, size_1 = 3, magnitude = TRUE)

# positive effect increasing in magnitude

##################################################
# Figure 3, Panel A
##################################################

# plot data
plot_data <- tibble(
  x = factor(
    x = c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5), 
    levels = c(1, 2, 3, 4, 5), 
    labels = c("27 to 15", "15 to 7", "7 to 5", "5 to 3", "3 to 1")
  ),
  y = c(
    calculate_cs(sim_0 = sim_baseline, size_0 = 27, size_1 = 15),
    calculate_cs(sim_0 = sim_judges, size_0 = 27, size_1 = 15),
    calculate_cs(sim_0 = sim_baseline, size_0 = 15, size_1 = 7),
    calculate_cs(sim_0 = sim_judges, size_0 = 15, size_1 = 7),
    calculate_cs(sim_0 = sim_baseline, size_0 = 7, size_1 = 5),
    calculate_cs(sim_0 = sim_judges, size_0 = 7, size_1 = 5),
    calculate_cs(sim_0 = sim_baseline, size_0 = 5, size_1 = 3),
    calculate_cs(sim_0 = sim_judges, size_0 = 5, size_1 = 3),
    calculate_cs(sim_0 = sim_baseline, size_0 = 3, size_1 = 1),
    calculate_cs(sim_0 = sim_judges, size_0 = 3, size_1 = 1)
  ),
  z = factor(
    x = c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1),
    levels = c(0, 1),
    labels = c("Low heterogeneity", "High heterogeneity")
  )
)

# make plot
panel_a <- ggplot(plot_data) + 
  geom_point(aes(x = x, y = y, color = z), size = 3) + 
  geom_segment(aes(x = x[1], y = y[1], xend = x[2], yend = y[2]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[3], y = y[3], xend = x[4], yend = y[4]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[5], y = y[5], xend = x[6], yend = y[6]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[7], y = y[7], xend = x[8], yend = y[8]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[9], y = y[9], xend = x[10], yend = y[10]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_point(aes(x = x, y = y, color = z), size = 3) + 
  scale_color_manual(values = c("gray70", "black"), name = NULL) + 
  theme_minimal() + 
  theme(legend.position = "bottom", legend.direction = "vertical") +
  titles_minimal(title = "Panel A", x = "Change in panel size", y = "Effect of change\non the average variance")

##################################################
# Figure 3, Panel B
##################################################

# plot data
plot_data <- tibble(
  x = factor(
    x = c(1, 2, 3, 4, 5), 
    levels = c(1, 2, 3, 4, 5), 
    labels = c("27 to 15", "15 to 7", "7 to 5", "5 to 3", "3 to 1")
  ),
  y = c(
    calculate_cs(sim_0 = sim_baseline, sim_1 = sim_judges, size_0 = 27, size_1 = 15),
    calculate_cs(sim_0 = sim_baseline, sim_1 = sim_judges, size_0 = 15, size_1 = 7),
    calculate_cs(sim_0 = sim_baseline, sim_1 = sim_judges, size_0 = 7, size_1 = 5),
    calculate_cs(sim_0 = sim_baseline, sim_1 = sim_judges, size_0 = 5, size_1 = 3),
    calculate_cs(sim_0 = sim_baseline, sim_1 = sim_judges, size_0 = 3, size_1 = 1)
  )
)

# make plot
panel_b <- ggplot(plot_data) + 
  geom_point(aes(x = x, y = y), color = "black", size = 3) + 
  geom_segment(aes(x = x[1], y = -Inf, xend = x[1], yend = y[1]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[2], y = -Inf, xend = x[2], yend = y[2]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[3], y = -Inf, xend = x[3], yend = y[3]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[4], y = -Inf, xend = x[4], yend = y[4]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[5], y = -Inf, xend = x[5], yend = y[5]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_point(aes(x = x, y = y), color = "black", size = 3) + 
  scale_y_continuous(limits = c(0, 0.018)) + 
  theme_minimal() + 
  titles_minimal(title = "Panel B", x = "Change in panel size", y = "Difference in the effect of change\non the average variance")

##################################################
# Figure 3
##################################################

# combine panels
figure_3 <- panel_a + panel_b + plot_layout(ncol = 2)

# save plot
ggsave(plot = figure_3, filename = "replication-materials/figure-3.pdf", device = "pdf", width = 10, height = 6, scale = 1.25)

##################################################
# judge selection (right)
##################################################

calculate_cs(sim_0 = sim_baseline, size_0 = 5, size_1 = 3) # cut-point mean at 0
calculate_cs(sim_0 = sim_case_facts, size_0 = 5, size_1 = 3) # cut-point mean at 1
calculate_cs(sim_0 = sim_baseline, sim_1 = sim_case_facts, size_0 = 5, size_1 = 3, magnitude = FALSE)
calculate_cs(sim_0 = sim_baseline, sim_1 = sim_case_facts, size_0 = 5, size_1 = 3, magnitude = TRUE)

# positive effect decreasing in magnitude 

##################################################
# Figure 4, Panel A
##################################################

# plot data
plot_data <- tibble(
  x = factor(
    x = c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5), 
    levels = c(1, 2, 3, 4, 5), 
    labels = c("27 to 15", "15 to 7", "7 to 5", "5 to 3", "3 to 1")
  ),
  y = c(
    calculate_cs(sim_0 = sim_baseline, size_0 = 27, size_1 = 15),
    calculate_cs(sim_0 = sim_case_facts, size_0 = 27, size_1 = 15),
    calculate_cs(sim_0 = sim_baseline, size_0 = 15, size_1 = 7),
    calculate_cs(sim_0 = sim_case_facts, size_0 = 15, size_1 = 7),
    calculate_cs(sim_0 = sim_baseline, size_0 = 7, size_1 = 5),
    calculate_cs(sim_0 = sim_case_facts, size_0 = 7, size_1 = 5),
    calculate_cs(sim_0 = sim_baseline, size_0 = 5, size_1 = 3),
    calculate_cs(sim_0 = sim_case_facts, size_0 = 5, size_1 = 3),
    calculate_cs(sim_0 = sim_baseline, size_0 = 3, size_1 = 1),
    calculate_cs(sim_0 = sim_case_facts, size_0 = 3, size_1 = 1)
  ),
  z = factor(
    x = c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1),
    levels = c(0, 1),
    labels = c("Neutral case facts", "Pro-plaintiff case facts")
  )
)

# make plot
panel_a <- ggplot(plot_data) + 
  geom_point(aes(x = x, y = y, color = z), size = 3) + 
  geom_segment(aes(x = x[1], y = y[1], xend = x[2], yend = y[2]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[3], y = y[3], xend = x[4], yend = y[4]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[5], y = y[5], xend = x[6], yend = y[6]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[7], y = y[7], xend = x[8], yend = y[8]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[9], y = y[9], xend = x[10], yend = y[10]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_point(aes(x = x, y = y, color = z), size = 3) + 
  scale_color_manual(values = c("gray70", "black"), name = NULL) + 
  theme_minimal() + 
  theme(legend.position = "bottom", legend.direction = "vertical") +
  titles_minimal(title = "Panel A", x = "Change in panel size", y = "Effect of change\non the average variance")

##################################################
# Figure 4, Panel B
##################################################

# plot data
plot_data <- tibble(
  x = factor(
    x = c(1, 2, 3, 4, 5), 
    levels = c(1, 2, 3, 4, 5), 
    labels = c("27 to 15", "15 to 7", "7 to 5", "5 to 3", "3 to 1")
  ),
  y = c(
    calculate_cs(sim_0 = sim_baseline, sim_1 = sim_case_facts, size_0 = 27, size_1 = 15),
    calculate_cs(sim_0 = sim_baseline, sim_1 = sim_case_facts, size_0 = 15, size_1 = 7),
    calculate_cs(sim_0 = sim_baseline, sim_1 = sim_case_facts, size_0 = 7, size_1 = 5),
    calculate_cs(sim_0 = sim_baseline, sim_1 = sim_case_facts, size_0 = 5, size_1 = 3),
    calculate_cs(sim_0 = sim_baseline, sim_1 = sim_case_facts, size_0 = 3, size_1 = 1)
  )
)

# make plot
panel_b <- ggplot(plot_data) + 
  geom_point(aes(x = x, y = y), color = "black", size = 3) + 
  geom_segment(aes(x = x[1], y = 0, xend = x[1], yend = y[1]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[2], y = 0, xend = x[2], yend = y[2]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[3], y = 0, xend = x[3], yend = y[3]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[4], y = 0, xend = x[4], yend = y[4]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_segment(aes(x = x[5], y = 0, xend = x[5], yend = y[5]), linetype = "dashed", size = 0.25, color = "gray80") +
  geom_hline(yintercept = 0, color = "black", size = 0.5) +
  geom_point(aes(x = x, y = y), color = "black", size = 3) + 
  theme_minimal() + 
  titles_minimal(title = "Panel B", x = "Change in panel size", y = "Difference in the effect of change\non the average variance")

##################################################
# Figure 4
##################################################

# combine panels
figure_4 <- panel_a + panel_b + plot_layout(ncol = 2)

# save plot
ggsave(plot = figure_4, filename = "replication-materials/figure-4.pdf", device = "pdf", width = 10, height = 6, scale = 1.25)

################################################################################
# end R script
################################################################################
